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We study the ground state phase diagram of the quantum spin-2 XXZ chain in the presence of on- 
site anisotropy using a matrix-product state based infinite system density-matrix-renormalization- 
group (iDMRG) algorithm. One of the interests in this system is in connecting the highly quantum 
mechanical spin-1 phase diagram with the classical S = oo phase diagram. Several of the recent 
advances within DMRG make it possible to perform a detailed analysis of the whole phase diagram. 
We consider different types of on-site anisotropies which allows us to establish the validity of the 
following statements: One, the spin-2 model can be tuned into a phase which is equivalent to the 
"topologically nontrivial" spin-1 Haldane phase. Two, the spin-2 Haldane phase at the isotropic 
Heisenberg point is adiabatically connected to the "trivial" large- D phase, with a continuous change 
of the Hamiltonian parameters. Furthermore, we study the spin-3 XXZ chain to help explain the 
development of the classical phase diagram. We present details on how to use the iDMRG method 
to map out the phase diagram and include an extensive discussion of the numerical methods. 



I. INTRODUCTION 

Quantum spin chains have proven to be extremely use- 
ful model systems for the study of strongly correlated 
quantum systems. In particular, many different types of 
phases and phase transitions can be understood by study- 
ing relatively simple Hamiltonians. A prime example is 
the SU(2) symmetric Heisenberg chain; it has gapless ex- 
citations for half integer spins, while the ground state of 
integer spin chains is protected by a gap in the energy 
spectrum. The existence of this gap was predicted by 
Haldane in the early eighties in Refs. Q] and [21 Follow- 
ing Haldane's predictions, Affleck, Kennedy, Lieb, and 
Tasaki (AKLT) presented model Hamiltonians for which 
the ground state can be obtained exactlyP^ The AKLT 
state was later found to exhibit unexpected properties, 
such as a nonlocal "string order" and edge states, which 
extend to other states within the same phasePEl The 
phase including the SU(2) symmetric point has conse- 
quently been referred to as the "Haldane phase" . 

These gapped ground states do not break any symme- 
try and thus cannot be characterized by any local or- 
der parameter. Instead, they are characterized by the 
projective representations of the symmetries present 
Physically, this means that the spin fractionalizes into 
two half-integer edge spins in the case of odd integer 
spin chains (one on each edge) , and into two integer edge 
spins in the even integer spin chains. There is a crucial 
difference between these two casesP^ In the odd integer 
case, the Haldane phase is a symmetry-protected topologi- 
cal phase (SPTP) as the half-integer edge spin cannot be 
removed unless the system undergoes a phase transition 
or all the relevant symmetries are explicitly broken. In 
contrast, in the even case the integer spins at the edge 
are not protected and thus the ground state can be adi- 
abatically turned into a trivial (product) state, without 



breaking any symmetries. This motivates the notation of 
two distinct phases, an odd-Haldane (OH) and an even- 
Haldane (EH) phase, which we adopt here. Seemingly, 
this would suggest the absence of a topological phase in 
the even integer spin chains. However, as proposed by 
OshikawaP all the Haldane phases corresponding to lower 
integer spin can in principle be realized in the presence 
of on-site anisotropy. In the spin-2 case, in particular, 
this implies the presence of the OH phase. It is one of 
the main goals of this work to verify if, and under what 
conditions, this scenario is realized. 

The study of such questions has been greatly ad- 
vanced by the development of matrix-product-states 
(MPS)^"2S1 and the reinterpretation of density-matrix- 
renormalization-group (DMRG) algorithms in terms of 
these states)2IH23] Infinite-system algorithms directly ob- 
tain the thermodynamic limit of infinite system size, 
without any finite-size corrections. This is especially cru- 
cial in the vicinity of critical phases, as will turn out to be 
the case in this work, with diverging correlation lengths. 
We adopt the infinit e-syst em DMRG, or iDMRG, algo- 
rithm in this paperPSHH Although this algorithm has 
been discussed priorly in the literature, we believe it is 
beneficial to give a detailed, pedagogical account of the 
algorithm and to compare and contrast it to the related 
infinite-system time-evolving-block-decimation (iTEBD) 
algorithmPSMlThis constitutes another main goal of this 
work. 

The remainder of the paper is organized as follows: 
We begin by presenting the model and reviewing rele- 
vant prior work in Sec. |Hj We present our spin-2 phase 
diagram and summarize its most important aspects in 
a summary of our main results in Sec. Ill In Sec. IV 



we formulate the iTEBD and iDMRG algorithms using 
a consistent notation. This allows us to highlight the 
similarities and differences of the algorithms. In Sec. [V] 
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we review the phases and phase transitions present in the 
model and discuss how they are most accurately observed 
with iDMRG. Details of our results, including data for 
an SPTP in a spin-2 chain and the adiabatic connection 
of the AFM Heisenberg point and points at large-D, are 
presented in the same section. We wrap up our investiga- 
tion by presenting the first, to our knowledge, numerical 
results for the spin-3 XXZ-chain, in Sec. |VI| 



II. MODEL 
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A. Definitions and prior studies 

The model Hamiltonian we are concerned with in this 
work is of the form 



H = H 



xxz 



D ■ 



(1) 



The first term describes the XXZ quantum spin chain 

^ {^n^n+i + Sn^n+i + AS*S* +1 ) , (2) 



H 



xxz 



where S", with a = x,y, z, is the a-component of the 
spin-5 operator at site n. A is the XXZ anisotropic inter- 
action parameter. As A is tuned from — oo to oo, the fol- 
lowing phases appear: ferromagnetic (FM), XY, Haldane 
(OH/EH, for integer spins only) and an ti-ferromagnetic 
(AFM), see Fig. [lj for S = ±, 1, 2.EEM1 In the following 
we will set the overall energy scale to J = 1. The transi- 
tion into the FM phase occurs at A = —1 for all S. The 
Heisenberg point A = 1 is the transition into the AFM 
phase for half integer S. For integer S it is within the 
Haldane phase, whose width in A decreases rapidly with 
increasing S consistent with its absence in the classical 
S — > oo phase diagram!^ This makes it computationally 
demanding to resolve the Haldane phases for large S. 

The term Hjj represents an on-site anisotropy, which 
for spin 5* has the general form 



2.S' 



(3) 



with D p constants. Depending on the values of the D p , 
this on-site anisotropy can favor all the different eigen- 
states \m a ) of a nonintcracting spin 5, and allows for 
realization of more phases than the XXZ model alone. 
The terms with odd power tend to favor the magneti- 
cally ordered AFM and FM phases, so we do not include 
them here. In the case of spin-2, there are then two terms 
with coupling constants D 2 and D4. A large positive D 2 
(or a large D4) favors a product state of the form 



\4>d) =® |0)„, 



(4) 



where |0) n is the eigenstate of with eigenvalue m z = 
0. This is because the higher m z eigenstates have been 



FIG. 1. (Color online) The phases present in a spin-S 1 XXZ- 
chain, for S = §,1,2. The Haldane phases appear only for 
integer spins and their width decreases rapidly with increasing 
S. At the same time, the phase transitions into the Haldane 
phases approach A = 1, where the direct XY-AFM phase 
transition occurs for all half-integer S. 



energetically projected out, and there remains only an 
effective spin-0 degree of freedom. In the limit D4 —> 
00 with D2 = — D4, only the states with m z — ±2 are 
projected out, and an effective spin-1 degree of freedom 
remains. This allows for the exploration of spin-1 phases, 
such as the OH phase, in the spin-2 chain. While nonzero 
D 2 has been considered before, see Ref . [3TH5SI we are not 
aware of any study that also includes the higher order D4 
term. 

We start by revisiting the case with D 2 as the only 
nonzero on-site anisotr opy. For spin-1, the phase dia- 
gram is well establisheoP^^ and is schematically plot- 
ted in Fig. [2j With increasing D2, there is a transition 
into a phase that contains the trivial product state Q 
as a ground state. This phase is often called the 'large- 
D' phase. However, in light of the recent classification 
scheme in terms of SPTP^Hni the large-D phase is ex- 
pected to be representative of the EH phase, so we will 
denote it as such (in fact, it can be thought of as the 
trivial S — Haldane phase). 

The situation for S > 2 is not as clear. Bosonization 
predicts that the spin-1 phase diagram is representative 
of all integer spina^ for small D2. Oshikawa predicted 
intermediate phases, corresponding to all the phases at 
the Heisenberg point for lower integer spins, between the 
Haldane phase and the large-D phaseP Early DMRG 
studie*2"23] did no t corroborate this picture but found 
for S — 2 a phase diagram already close to the semi- 
classical limit S — >• 00 (inset Fig. [2]). The Haldane phase 
surrounding the Heisenberg point appeared to be sepa- 
rated from the trivial large-D phase by t he XY phase. 
Later level spectroscopy (LS) studies^^U suggested in- 
stead that the Haldane phase is connected with the large- 
D phase; furthermore, these studies concluded in favor 
of the presence of a tiny OH phase in the D4 = plane 
around (A = 2.4, D 2 = 2). 

The two main questions, that have been hard to verify 



FIG. 2. A sketch of the spin-1 phase diagram from Ref. [41] 
(with our labeling of the phases). Inset: Phase diagram in 
the semiclassical limit S — > oo. 



for the spin-2 case, and which we address here are: i) 
Is the Heisenberg point indeed adiabatically connected 
to the trivial large- D phase? ii) Is there an OH phase 
in the phase diagram, particularly in the experimentally 
relevant case of small on-site anisotropy? To answer these 
questions, we will find it useful to introduce nonzero -D4, 
as this will allow us to easily access the OH phase and 
determine its extent. 



III. SUMMARY OF THE MAIN RESULTS 

In this section we summarize our main S — 2 results by 
discussing the phase diagrams obtained in our iDRMG 
simulations and shown in Fig. [3] and Fig. [4j deferring 
the details of how they were obtained to Sec. [V] The 
phase diagram for the relevant upper right corner A > 
0, D 2 > at L> 4 = is shown in Fig. [3| There is no 
direct transition from the XY phase into the AFM phase. 
Instead, these phases are separated everywhere by the 
EH phase which is continuously present to large D 2 . This 
answers the first of our main questions; the Heisenberg 
point and the large- D phase are continuously connected. 

The extension of the XY phase is consistent with the 
recent LS study of Ref |551 but covers a much smaller 
range of A than obtained in the earlier DMRG study 
of Ref. 1331 However, our numerical results suggest the 
complete absence of an OH phase in this plane, in con- 
trast to Refs. 1551 and 1551 where it is observed next to the 
XY phase. We believe that the reason for these different 
conclusions arises from the proximity of an EH <-> OH 
phase transition at small D4, as discussed below, which 
makes the relevant parameter region appear critical in 
finite chains. 

The -D4 anisotropy is an important parameter in estab- 



FIG. 3. (Color online) The spin-2 phase diagram at D4 = 
for A > and D2 > as obtained with iDMRG. Examples 
of data used to obtain the phase boundaries are presented in 
Sec.jV] In particular, data corresponding to the vertical line is 
given in Figs. |12| and |13| and the tilted line in Figs. |13| and |14| 



lishing the above result. In fact, since large D 4 = —D 2 
effectively realizes a spin-1 chain, one expects the OH 
phase to appear at the Heisenberg point along the lines 
D 2 = -D4 + Z)f =1 witn -°2 =1 tne values at which the 
OH phase appears in the S — 1 phase diagram. This is 
indeed what we observe, as shown in Fig. [4] The fact that 
the OH phase is easily realized when introducing the ad- 
ditional anisotropy parameter D4 is one of the main new 
result of our work. The question of whether it is obtained 
in the -D4 = plane is then simply the question of how 
large the extent of the OH phase is. While we believe 
that our results demonstrate that it does not touch the 
-D4 = plane, we cannot strictly rule out that it does, as 
discussed in more detail in Sec. [V] 

IV. NUMERICAL TECHNIQUES 

In this section we outline the numerical method used to 
obtain our results. Though most of the details have been 
discussed in the literature elsewhere, we find it useful to 
discuss our implementation of the algorithm. The infinite 
time evolving block decimation (iTEBD^pH anc i the infi- 
nite density matrix renormalization group (iDMRG) 24 
algorithms are both based on the infinite matrix-product 
state (iMPS) representation. 2 - As we explain shortly, 
MPSs can efficiently represent many-body wave functions 
where the accuracy is controlled by the bond dimension 
X (the error decreases rapidly with increasing \)- Using 
methods which work with infinite systems has a number 
of advantages: no extrapolation to the thermodynamic 
limit is needed; there are no edge modes which can com- 
plicate the convergence of the algorithm; and, as shown 
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FIG. 4. (Color online) The spin-2 phase diagram at A = 1.0 
as obtained with iDMRG. The OH phase is easily reached 
by including the D4 on-site anisotropy term. The dots in- 
dicate the values at which the simulations were performed. 
To resolve the narrow EH phase, simulations on a finer grid 
were performed between the AFM and XY phases (points not 
shown) . 



later on in this section, finite entanglement scaling can 
be used to extract quantities such as the central charge. 
We begin by reviewing some details of this infinite sys- 
tem representation focusing on translationally invariant 
systems, and then contrast and compare the two numer- 
ical methods using a consistent notation. We do not aim 
to provide a complete discussion of the techniques but 
rather a clear and compact introduction to the methods 
used. 

The concept of entanglement is central to the MPS rep- 
resentation and the algorithms based on it. The so-called 
entanglement spectrum^ is obtained from the Schmidt 
decomposition (singular value decomposition): Given a 
bipartition T~L = Hl ®T~Lr of the Hilbert space (below 
%l and Hr, respectively, represent the states on sites to 
the left and right of a bond), any state |>P) g W can be 
decomposed as 

I*) =^A Q |a) L «> \a) R , \a) R/L €H R/L - (5) 

a 

The Schmidt coefficients (singular values) A Q can al- 
ways be chosen positive, the states {|o;)l} and 
form orthonormal sets in Hl and Hr respectively, i.e., 
(a\j3) L = (a\/3) R — 5 a p, and by normalization J2 a ^a = 
(^f\^) = 1. The Schmidt decomposition is related to 
the reduced density matrix for one half of the system, 
p R = Tr-Hj- (|^)(^|). In particular, the Schmidt states 
are the eigenstates of p R and the Schmidt coefficients are 
the square roots of the corresponding eigenvalues, i.e., 



p R = J2 a ^a\ a )( a \R ( an d analogously for p L ). This di- 
rectly gives the entanglement entropy through 

(6) 



S £ = -]TA>gA2. 



Finally, the entanglement spectrum {e Q } is related to the 
spectrum {A^} of the bipartition by A^, = exp(— e a ) for 
each a. 



A. Matrix-product states 

A general quantum state |\&) on a chain with N sites 
can be written in the following MPS formliSHini 

1*)= Y, A^AW*...A™»\j 1 ,...,j N ). (7) 

Here, A^ Jn is a Xn-i x Xn matrix and \j n ) with j n = 
1, . . . , d is a basis of local states at site n. We call the 
indices of the matrices "bond" indices. The matrices at 
the boundary, i.e., n = 1 and n — N, are vectors, that is 
Xo = Xn = 1) such that the matrix product in ([7| pro- 
duces a scalar coefficient. The superscript [n] denotes the 
fact that for a generic state, each site is represented by 
a different set of matrices. Ground states of one dimen- 
sional gapp ed sy stems can be efficiently approximated 
by an MPSj^El in the sense that the value of the %'s 
needed to approximate the ground state wave function 
to an arbitrary precision is finite as N —¥ 00. The phys- 
ical insight that allows us to make this statem ent is the 
area law, which holds for this class of systems? 45 * 46 ! De- 
tails on how the accuracy of the representations depends 
on x can be found in Ref. HU 

Canonical form. Without a loss of generality, we 
write the matrices A J as a product of Xj-i x Xj complex 
matrices T 3 and positive, real, square diagonal matrices 
A, 



I*) 



E 

jl,—,jN 



rMjiAWrPi^AP] • • • A^-^rM™ 

X |ji,...,jjv) 



(8) 



which is pictorially illustrated in Figs. [5^ a) and[5]jb). A 
rank-n tensor is represented by a symbol with n protrud- 
ing lines. (For example, T, a rank-3 tensor, has three 
indices and is represented by a triangle with three lines 
protruding from it.) Connecting the lines among tensors 
symbolizes a tensor contraction, i.e., summing over the 
relevant indices. In the following we will motivate the 
choice (J||| for the MPS form. 

Equation ^ allows for many possible representation 
of the same wave function, as we can insert a resolution of 
the identity 1 = XX^ 1 into any bond. This freedom can 
be used to define a 'canonical form' of the MPS, following 
Ref. [551 and 27J Any bond n defines a bipartition of the 
system into sites L — {1, . . . , n} and R = {n+ 1, . . . , N} 
to the left and right of the bond. From the form of the 
MPS, we can define a set of Xn wave functions la)^^ to 
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FIG. 5. Diagrammatic representation of (a) the tensors 
F and A. The horizontal lines represent the bond indices 
a G {1, ...,x} and the vertical lines the physical indices 
j € {1, . . . , d}. (b) An MPS formed by the tensors T and 
A. Connected lines between tensors (or within a single ten- 
sor) denote summation over the corresponding indices, (c) 
Definition of the right and left (Schmidt) basis states with 
respect to a partition on a bond with index a. (d) Condition 
for the MPS to be in the canonical form. The transfer matrix 
T of Eq. ( |10[ | has been shaded. The upside-down triangles 
are the complex conjugate of the V tensors, (e) If the state is 
in canonical form, then the dominant right eigenvector of T 
is the 'identity matrix' with eigenvalue equal to 1. A similar 
condition applies for the left transfer matrix T. (f) The cor- 
relation function defined in Eq. (131. The squares correspond 
to the operators F(0) and G(r). 



the left/right of the bond [see Fig. [5jc)] such that state 
takes the form 



(9) 



The MPS representation {rW , A^ , . . . , T^} is in canon- 
ical form if: For every bond, the set of wave functions 
\oi)^ R along with A^ form a Schmidt decomposition of 

. In other words we must have (a\a)^ — <5 QQ and 

(a\a)^ — 5 &a , along with ^(A Q *') 2 = 1 on every bond. 
For finite systems, a generic MPS can be transformed 
into canonical form by successively orthogonalizing the 
bonds starting from either the left or right end of the 



chainPS 

Infinite matrix product states. In this paper we are 
most interested in infinite chains, N — > oo. In this case, 
translational invariance is restored and the set of matri- 
ces on any given site becomes the same, that is rWJ = 
and A[™1 = A for all integers n. To check if the iMPS is in 
canonical form, we need to compute the overlaps (a\a)n, 
which would appear to require an infinite tensor contrac- 
tion. But, we can use the translation invariance to pro- 
ceed inductively. For infinite MPS, we can conveniently 
express the orthogonality condition (i.e., canonical form) 
in terms of a transfer matrix T [illustrated in Fig. ^d)] 
defined as 



T. 



aa;P0 



A/3A^, 



(10) 



where '*' denotes complex conjugation (and is pictorially 
represented by an upside-down triangle). The transfer 
matrix T relates the overlaps defined on bond n with 
overlaps defined on bond n + 1. Given that the right 



basis states 
states |a)jj 



IR 



i+l] 



on bond n + 1 are orthonormal, the 

on bond n will also be orthonormal if T has 
a right eigenvector Sa§(= 1) with eigenvalue 77 = 1, as 
illustrated in Fig. [SJe) . For the left set of states we define 
an analogous transfer matrix T, 



T 



3 



A„A, 



1 a/3 



(T J -) 



(11) 



which must have a left eigenvector 5 a& with rj = 1. These 
eigenvector criteria are clearly necessary conditions for 
all bonds to be canonical; in fact, assuming in addition 
that 77 = 1 is the dominant eigenvalue, they are sufficient. 
An algorithm to explicitly transform an arbitrary infinite 
MPS to the canonical form is given in Ref. |3H1 

If the infinite MPS is not translational invariant with 
respect to a one-site unit cell, all the above can be simply 
generalized by considering a unit-cell of L sites which 
repeats itself, e.g., in the case of L = 2 the tensors are 
given by 



p[2n] 
p[2n+l] 



= r A , 
= r B , 



A [2n]j =A/lj 

A.pn+1] = A B 



(12) 



for n £ Z. Reviews of MPS's as well as the canonical 
form can be found in Refs. [2H1 HI and 1391 

Calculations of observables from an iMPS. If the 
MPS is given in canonical form, we can use the orthog- 
onality of the Schmidt states to evaluate local expecta- 
tion values by contracting the tensors locallypS Corre- 
lation functions can be obtain using the transfer matrix 
Eq. ([10]). For this we evaluate (F(O)G(r)) of an iMPS. 
Let r > 0, then 

(F(O)G(r)) =T L (F)T r - 1 T R (G), 

T L (F) aa = ^A2(F 7a )*F«F 7a A a A a , 
7 

T *XG0^£A 7 (ry*G^ r (13) 
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Ti ir are the "stubs" which measures F and G locally, in 
between which we put r — 1 copies of the transfer matrix 
T, see Fig. [5jf ) for a pictorial representation. Local ob- 
servables (F[0)) can be obtained from the same expres- 
sion, replacing Tr(G) and T° with identity operators. 
By generalizing the transfer matrix to include on-site op- 
erators, non-local order parameters can be obtained with 
the same approach 



r A A A r B A B r A A A r B A B r A A A 

■■^j^o-^-o-^-o-^-o-^- ~ 



u 



u 



u 



u 



u 



T, 



R 



E 



(r 



AaA 



(14) 



For example, calculating the correlation function with 



F = G = S; 

Eq. 



(44) 



and R = e i7rS ~ , we obtain the "string order" 



The resulting correlation functions generically take the 
form of a sum of exponentials, with the slowest decaying 
exponential determined by the second largest (in terms 
of absolute value) eigenvalue ti of the transfer matrix. 
We define the correlation length of the MPS as 



i 



log |e 2 | 



(15) 



which is readily obtained using a sparse algorithm to 
find the eigenvalues of the transfer matrix. A degen- 
erate largest eigenvalue indicates that the state is in a 
'cat state,' i.e., in a superposition of different superselec- 
tion sectors, which can occur when there is spontaneous 
symmetry breaking!^ 

In systems with a conserved quantum number (e.g., the 
total S z ), one can calculate the correlation length for op- 
erators (F, G) corresponding to different sectors from the 
corresponding eigenvalues of the transfer matrix. In this 
paper we denote the two correlation lengths correspond- 
ing to operators which change the quantum numbers by 
S z = (e.g. (S Z S Z )) as £ and S z = ±1 (e.g. (S+S - )) as 
£i. The correlation length £ is given by the largest one, 
i.e., £ = max(£ ,£i, • ■•)• 



B. Infinite Time Evolving Block Decimation 
(iTEBD) 

In the iTEBD algorithm, we are interested in evaluat- 
ing the time evolution of a quantum state: 



(16) 



The time evolution operator U can either be U(t) = 
exp(-iHt) yielding a real time evolution, or an imagi- 
nary time evolution U(t) — exp(— Ht). The latter is 
used to find ground states of the Hamiltonian H through 
the relation 



\1> G ) = lim e- TH |Vo> 



(17) 



To achieve this, one makes use of the Trotter-Suzuki de- 
composition, which approximates the exponent of a sum 



FIG. 6. In iTEBD each time step St of a time evolution is 
approximated using a Trotter-Suzuki decomposition, i.e., the 
time evolution operator is expressed as a product of unitary 
two-site operators. 



of operators, with a product of exponents of the same 
operators. For example, the first order expansion reads 



e (A+B)S = e AS e BS 



0(S 2 ). 



(18) 



Here A and B are operators, and 5 is a small parameter. 
The second order expansion similarly reads 



e AS/2 £ BS e AS/2 



0(S 3 



(19) 



To make use of these expressions, we assume that the 
Hamiltonian is a sum of two-site operators of the form 
H = ^2 n /j[ n ' rl + 1 ] and decompose it as a sum 



H — H m 

= £' 

n odd 



He, 



1 even 
[n,n+l] i 



n even 



(20) 



Each term i/ 0( jd and H cvcn consists of a sum of commut- 
ing operators. 

We now divide the time into small time slices St <£L 1 
(the relevant time scale is in fact the inverse gap) and 
consider a time evolution operator U(St). Using, as an 
example, the first order decomposition (18), the operator 



U(8t) can be expanded into products of two-site unitary 
operators 



U(St) 



where 



odd 



Yl U [n < n+1] {5t) 



U [n > n+1] {5t) = 



(21) 



(22) 



This decomposition of the time evolution operator is 
shown pictorially in Fig. [6] One notices that even if the 
underlying system has a translation invariance of one site, 
the decomposition breaks this temporarily into a two site 
translation symmetry. Therefore, one needs to keep at 
least two sets of matrices T A ,A A and T B ,A B . The suc- 
cessive application of these two-site unitary operators to 
an MPS is the main part of the algorithm. 
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(25) 
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FIG. 7. The iTEBD update scheme for a unitary two-site 
transformation of a two-site unit cell MPS in canonical form 
(see text for details). 



Local unitary updates of an MPS. One of the advan- 
tages of the MPS representation is that local transforma- 
tions can be performed efficiently. Moreover, the canon- 
ical form discussed above is preserved if the transforma- 
tions are unitary.^ 

A one-site unitary U simply transforms the tensors F 
of the MPS 



^ = E 



TJ3 yj 



(23) 



If we consider an infinite, translational invariant MPS, 
this transformations implies the application of the uni- 
tary to all equivalent sites simultaneously. In such case 
the entanglement of the wave-function is not affected and 
thus the values of A do not change. 

The update procedure for a two-site unitary transfor- 
mation acting on two neighboring sites is shown in Fig. [7] 
We focus on an update of an AB bond between two neigh- 
boring sites n and n + 1 for an MPS with a unit cell of 
size N = 2. The inequivalent BA bonds are updated 
similarly by simply exchanging A and B, The general- 
ization to an L-site unit cell is straightforward. We first 
find the wave function in the basis spanned by the left 
Schmidt states on bond n—l:n, the l-site Hilbert space 
of sites n and n + 1, and the right Schmidt states on 
bond n + 1 : n + 2, which together form an orthonormal 
basis {\a n -x) L , \j n ), \k n +i), |7™+i)flJ- Calling the wave 
function coefficients 9, the state is expressed as 

IV>> = E ei fe 7 K-i) i |in)|fc„ +1 )|7n+i)i i - (24) 

Using the definitions of \a) L / R shown in Fig. [5](b) , 9 is 



Writing the wave function in this basis is useful because 
it is easy to apply the two-site unitary in step (ii) of the 
algorithm: 



a 7 / j i'k f cry ' 
j'k< 



(26) 



Next we have to extract the new tensors f A , T B and A A 
from the transformed tensor 9 in a manner that preserves 
the canonical form. We first 'reshape' the tensor 9 by 
combining indices to obtain a d\ x d\ dimensional matrix 
9 JC[; fe 7 . Because the basis |a n _i) i\jn) is orthonormal, as 
for the right, it is natural to decompose the matrix using 
the singular value decomposition (SVD) in step (iii) into 
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(27) 



P 



where X, Y are isometries and D is a diagonal matrix. 
The isometry X relates the new Schmidt states |/3„)l 
to the combined bases \a n -i) L\j n ) ■ Analogously, the 
Schmidt states for the right site are obtained from the 
matrix Y. Thus the diagonal matrix D contains pre- 
cisely the Schmidt values of the transformed state, i.e., 
A A = D. The new tensors T A ,T B can be extracted di- 
rectly from the matrices X, Y using the old matrices A B 
and the definition of 9 in Eq. ( 25 ) . In particular we 
obtain the new tensors in step (iv) by 



f $ = (A B )«%-^ 

yB, 3 
L Pl 



(28) 
(29) 



After the update, the new MPS is still in the canonical 
form. Note that as in the one-site update, if we apply 
the algorithm to an MPS, the update is performed si- 
multaneously to all matrices at equivalent bonds. Thus 
the iTEBD algorithms exploits the translational invari- 
ance of the systems by effectively performing an infinite 
number of parallel updates at each step. 

The entanglement at the bond n, n + 1 has, in the up- 
date, changed and the bond dimension increased to d\- 
Thus the amount of information in the wave function 
grows exponentially if we successively apply unitaries to 
the state. To overcome this problem, we perform an ap- 
proximation by fixing the maximal number of Schmidt 
terms to After each step, only the x most important 
states are kept, i.e., if we order the Schmidt states ac- 
cording to their size we simply truncate the range of the 
index /3 in (27 1 to be 1 . . . \. This approximation lim- 



its the dimension of the MPS and the tensors T have at 
most a dimension of d x x x X- Given that the trun- 
cated weight is small, the normalization conditions for 
the canonical form will be fulfilled to a good approxi- 
mation. In order to keep the wave function normalized, 
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FIG. 8. (a) An operator O acting on an entire chain expressed 
as a matrix product operator, (b) A matrix product operator 
acting on a matrix product state 0\%p). (c) The expectation 
value {if)\0\ip} expressed in an MPO form. 



one should divide by the norm after the truncation, i.e., 
divide by J\f = 



IP) 13 



j, a, 7 I " a 7 | ' 

If we perform an imaginary time evolution of the state, 
the operator U is not unitary and thus it does not con- 
serve the canonical form. It turns out, however, that 
the successive Schmidt decompositions assure a good ap- 
proximation as long as the time steps are chosen small 
enoughP^ One way to obtain very accurate results is to 
decrease the size of the time steps successively™ 

The simulation cost of this algorithm scales as d 3 x 3 
and the most time consuming part of the algorithm is 
the SVD in step (hi). If the Hamiltonian has symme- 
tries, we can considerably accelerate this step by explic- 
itly conserving the resulting constants of motion. The 
anisotropic spin model we study has for example a global 
U(l) symmetry and conserves the total magnetization. 
Thus the matrix Qi a -j-y has a block-diagonal form and 
the SVD can be performed in each block individually, 
yielding a considerable speed up. See Refs. l5TH53"l for the 
details of the implementation of symmetries into the algo- 
rithm. Numerically, the algorithm can become unstable 
when the values of A become very small since the matrix 
has to be inverted in order to extract the new tensors in 
step (iv) of the algorithm. This problem can be avoided 
applying a slightly modified version of this algorithm as 
introduced by Hastings in Ref. I5H 



C. Matrix-Product Operators 

The iDMRG algorithm explained in the next section 
relies on expressing the Hamiltonian of the system in 



terms of matrix product operator (MPO). An MPO is 
a natural generalization of an MPS to the space of oper- 
ators. An operator in an MPO form, acting on a chain 
with L sites, is given by 



O 



ji 



wicftM [lb ' lj; M [2b2 ^ • ■ • M[ L ]^'W rig ht (30) 

x\jl,---,jL){ji,-.-,j' L \ , 



,3l 



where M 3nJ ^' are DxD matrices, and \ j n ), \ j' n ) represent 
local states at site n, as before. At the boundaries we 
initiate and terminate the MPO by the vectors v\ e ft and 

Wright ■ 

A pictorial representation of an MPO is given in 
Fig. ]8^a). The notation is very similar to the one for 
an MPS: the horizontal line corresponds to the indices 
of the virtual dimension and the vertical lines represent 
the physical states \j n ) (bottom) and (j' n \ top. The ad- 
vantage of the MPO is that it can be applied efficiently 
to a matrix product state as shown in Fig. [8jb). All 
local Hamiltonians with only short range interactions 
can be represented using an MPO of a small dimen- 
sion D. Let us consider, for example, the MPO of the 
anisotropic Heisenberg model (J2| in the presence of an 
on-site anisotropy. Expressed as a tensor product, the 
Hamiltonian takes the following form: 



+ S y <£> S v ® 1 <8> • • • ® 
+ AS Z ® S z <g> 1 ® • • • < 
+ [D 2 (S Z ) 2 + D 4 (S z f 



1 + 1 
1 + 1 
g>l + ... 

® 1 ® 1 



S x - 



)1 + .. 
1 + .. 



1 + 



The corresponding exact MPO has a dimension D 
and is given by 



(31) 
= 5 



MW = 



1 

s :r 
s y 

AS* 2 



0\ 

' 







(32) 



with 

Vleft 



\D 2 (S 2; ) 2 +D4(S i ) 4 S x S v S z 1/ 

(0, 0, 0, 0, 1) , ukght = (1, 0, 0, 0, 0) 



T 

(33) 



By multiplying the matrices (and taking tensor products 
of the operators), one can easily see that the product 
of the matrices does in fact yield the Hamiltonian (31 1. 



Further details of the MPO form of operators can be 
found in Refs. |23] and M 



D. Infinite Density Matrix Renormalization Group 
(iDMRG) 

We now discuss the infinite Density Matrix Renor- 
malization Group (iDMRG) algorithm. Unlike iTEBD, 
the iDMRG is a variational approach to optimizing the 
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(r A )*A A (r B )*A B (r B f^ 





A A (r A )* A A (r A )* a 3 (r B )* 



FIG. 9. Pictorial representation of a contraction of the left 
and right environments. The boundaries are initiated by the 
tensors R^ iBl ,a = <5 a ,aWri g ht;a and L%,^,a = &*,aifiaftja. 



MPS, but the algorithms have many steps in common. 
One advantage of the iDMRG is that it does not rely on 
a Trotter-Suzuki decomposition of the Hamiltonian and 
thus applies to systems with longer range interactions. 
We assume only that the Hamiltonian has been written 
as an MPO. Secondly, the convergence of the iDMRG 
method to the ground state is in practice much faster. 
This is particularly the case if the gap above the ground 
state is small and the correlation length is long. 

The schematic idea for the iDMRG algorithm is as fol- 
lows (see Fig. 10 1. Like in iTEBD, the state at each 
step is represented by an MPS. We variationally opti- 
mize pairs of neighboring sites to minimize the ground 
state energy (ip\H\ip), while keeping the rest of the chain 
fixed. To do so, at each step we represent the initial wave 
function using the two site tensor 0jpL (as previously 
defined in the iTEBD section), project the Hamiltonian 
into the space spanned by the basis set \ajkj3), and use 
an iterative algorithm (e.g. Lanczos) to lower the energy. 
Repeating this step for each pair, the wave function con- 
verges to the ground state. For simplicity, only the details 
of the algorithm with a unit cell of two sites, A and B, 
will be described below. 

Two-site update algorithm. We start by describing 
the update of an AB bond between two neighboring sites 
n and n + 1 (the update on a B A bond can be performed 
analogously by exchanging the role of A and B), and re- 
turn later to the initialization procedure. Step (i) is iden- 
tical to the first step in the iTEBD method; we contract 
the tensors for two neighboring sites to obtain the initial 
wave function OjpL. The orthonormal basis \aj/3k) spans 
the variational space \ip) — QijL\ajf3k) of the update, in 
which we must minimize the energy E = {tp\H\%p) in or- 
der to determine the optimal 0. Because H is written 
as an infinite MPO, it appears at first that to evaluate 
the energy we will have to contract an infinite number 
of tensors starting from left and right infinity, as illus- 
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FIG. 10. A pictorial representation of an iDRMG iteration 
step update. See text for details. 



trated in Fig. J8jc). For the sake of induction, however, 
suppose we have already done this contraction on the left 
through bond n — 1 : n, and on the right through bond 
n + 1 : n + 2. As illustrated in Fig. |9j the result of these 
contractions can be summarized in two three leg tensors 
we call the left and right "environments." The left envi- 
ronment L a&J a has three indices: the MPO index a, and 
the indices a, a corresponding to the bond indices of \tp) 
and Likewise, on the right we have R-yy, c - Each bond 
of the system has a similarly defined environment; for a 
unit cell of two, we have in total {L A ,L B },{R A ,R B }. 
These environments are nothing other than the MPO for 
the Hamiltonian projected into the space of left and right 
Schmidt states about each bond. 

With the environment in hand, we can project the 
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Hamiltonian into the orthonormal basis \ajjk); to min- 
imize the energy of we find the ground state of the 
X 2 cP x x 2 d 2 "Hamiltonian": 



H. 



^aa,a 1V1 ab IVI bc U 



k,k U B 

77 j c ' 



(34) 



To find this ground state, we use an iterative procedure 
such as Lanczos or Jacobi-Davidson at a cost of x 3 Dd 2 
per multiplication, as illustrated in step (ii), and obtain 
an improved guess for the wave function O and energy 
Eq. By using the initial wave function O as the start- 
ing vector for the minimization procedure, convergence 
is typically reached with only a couple of steps. This can 
be compared to the iTEBD optimization where we ob- 
tain a new wave-function after applying the imaginary 
time-evolution operator. As with iTEBD, the bond di- 
mension grows as x ~ * d\, which we must truncate using 
SVD, shown in step (iii). It is important that the left 
and right Schmidt basis about any bond remain orthogo- 
nal, because we assume \aj(3k) is an orthogonal basis at 
each step. Assuming this was the case on bonds of type 
B, the isometry properties of the SVD matrices X and 
Y imply that the orthogonality condition holds for the 
updated Schmidt states defined about the central bond 
A, and hence will remain so throughout the simulation. 
At this point, we have improved guesses for the matrices 
f A / B ,A A in step (iv). 

The last step is to update the environment. At a min- 
imum, we must update the environments on the bond 
which we just optimized by simply multiplying new ten- 
sors to the left and right as shown in Fig. 10 V): 



= tB iflfi 



n 0f3,b — n il,a L /3 7 fe A 7 1V1 ab 1 7 /3fe 7V 7 ' 



(35) 
(36) 



This concludes the update on bond AB and we move over 
by one site, exchanging the roles of A and B, and repeat 
until convergence is reached. 

Initializing the environment. We now return to the 
problem of initializing the algorithm. The initial MPS 
can be arbitrary (though it should be in canonical form) . 
A fine choice is a x = 1 tensor product state which either 
preserves or breaks the symmetries as desired. To form 
the initial environment, we suppose when computing the 
left /right environment that H is zero to the left/right of 
the bond, which is captured by tensors of the form 



nH _ r _ ^ 



T m - — A - n 
^a,a,a u a,a u \cft;a> 



(37) 
(38) 



where the «i e ft /right are as m Eq. (30). Referring to Eq 



(33) as an example, recall that w r i g ht specifies the MPO 
index such that no further operators will be inserted to 
its right; likewise, v\ c ft indicates no operators have been 
inserted to its left. Because all terms in the Hamiltonian 
then act as the identity to the left /right of the bond, the 



orthogonality of the Schmidt vectors implies that pro- 
jecting the identity operator into the left/right Schmidt 
basis trivially gives 5 a ^ & . When symmetry breaking is ex- 
pected it is helpful to further initialize the environments 
by repeatedly performing the iDMRG update without 
performing the Lanczos optimization, which builds up 
environments using the initial symmetry broken MPS. 

Ground state energy from iDMRG. One subtlety of 
the above prescription lies in the interpretation of the 
energy Eqs obtained during the diagonalization step. Is 
it the (infinite) energy of the infinite system? Using the 
initialization procedure just outlined, the Lanczos energy 
Egs after the first step is the energy of the two-site prob- 
lem. While we motivated the environments as represent- 
ing infinite half chains, it is more accurate to assign them 
a length of after the initialization procedure, and at 
each optimization step the length of the left /right envi- 
ronment about the central bond increases because a site 
has been appended. Keeping track of the length £r/l of 
each environment (for a unit cell of two, each grows on 
alternate steps), we see that the energy Eqs corresponds 
to a system of size I = £l + 2 + £r. By monitoring the 
change in Eqq with increased i, we can extract the en- 
ergy per site. This is convenient for problems in which 
there is no few-site Hamiltonian with which to evaluate 
the energy. 

As for the iTEBD algorithm, we can achieve a consid- 
erable speed-up by using the symmetries of the Hamilto- 
nian, which requires assigning quantum numbers to the 
tensors of the MPO in addition to the MPS. 



E. Finite entanglement scaling 

An advantage of the infinite system methods intro- 
duced above is that no artifacts from the boundary ap- 
pear. On the other hand, finite size effects can be very 
useful for performing a scaling analysis. In this section 
we show that critical properties of the system can be ex- 
tracted by performing a "finite-entanglement scaling" in 
the infinite systems. This means, one can perform sim- 
ulations with different bond-dimensions x at a critical 
point and use the induced finite correlation length £ x as 
a scaling variable analogous to a finite system size. 

To motivate this notion, consider the entanglement en- 
tropy Se, which for an infinite system diverges logarith- 
mically as a function of the correlation length as critical- 
ity is approached.^ In an MPS, however, Se is bounded 
by Se < log Xi and an infinite x is needed to accurately 
represent critical states. Clearly we cannot perform sim- 
ulations with an infinite x^ raising the question: what 
happens if we nevertheless optimize a finite dimensional 
MPS for a critical system? This question has been ad- 
dressed by a series of papers EMU it turns out that simu- 
lating critical systems using finite x cuts off long distance 
correlations a finite length £ x . If we define the correlation 
length of the MPS £ to be the length obtained from the 
second largest eigenvalue of the transfer matrix, as define 
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in Eq. (15), then at criticality the correlation length of 



the MPS scales as 



where k 



(39) 



Because x introduces a length 



scale in a universal way, we can define the 'finite entan- 
glement length' by £ x = C\ K , where C is independent of 
X, and extract various quantities of interest using a finite 
£ x scaling analysis, or "finite entanglement scaling." In 
an infinite system at criticality, the scaling relations are 
generally obtained from the analogous scaling relations 
in a finite size system by replacing the finite length L by 
£ x . For example, for a critical point with central charge 
c, the entanglement entropy Se between two halves of a 
finite system of length L scales as Se = § log(L/a) + s , 
with a the lattice spacing and sq a non-universal con- 
stant. If we instead measure Se for an infinite system, 
but with finite x, we can substitute L — > £ x , 



Se = glogfo/a) + 4 



(40) 



The additive constant is again non-universal, and unre- 
lated to So. One should note, though, that while £ x and 
L have the same scaling dimension (i.e., that of length), 
the actual scaling functions are not guaranteed to be the 
same (for example, s 7^ s 'o m the above). 

Another useful quantity at criticality for systems with 
a U(l) symmetry is the 'stiffness,' here parameterized 
as the Luttinger parameter For Hamiltonians that 
conserve the total magnetization, it can be obtained from 
the scaling of bipartite spin fluctuations of a half chain 



F = ((Si) 2 ) ~ (Si)'' 



(41) 



where 5£ is the z-component of the total spin to the left 
of a cut (for example, the total mag netiza tion of the sites 
i < 0). The spin fluctuations satisfy^SEDl 



K 
2^2 



log(£ x /a) + const, 



(42) 



allowing us to extract K by measuring the scaling of F 
with increased x- 

The above discussion holds at criticality. Next we dis- 
cuss the situation in the vicinity of a critical point, where 
the physical correlation length £ p hys is finite but much 
larger than the correlation length induced by finite Xi 
i.e. £ x <C £ P hys- For our purposes, we can define £ p hys by 
the MPS correlation length £ in the limit x ~ * 00 1 where 
the MPS represents the true ground state. In this regime 
the MPS is still cutoff by £ x , rather than the true cor- 
relation length £ph ys> so the finite-entanglement scaling 
relations Eqs. (40 1 and (42) can still be used to obtain 



quantities like c and K. We will refer to this parameter 
range as the finite entanglement scaling region and it is 
shown as the red area in the schematic plot in Fig. [TT] 

Further away from criticality, as £ p hys we can 

fully converge the MPS and the state 'knows' it is not at 




FIG. 11. (Color online) A schematic plot illustrating the idea 
of finite-entanglement scaling. The black solid line shows the 
exact correlation length £ p h ys of the Hamiltonian, which di- 
verges at the critical point. The dotted lines show the corre- 
lation lengths £ of the optimized iMPS at a finite %2 > Xi- 
The horizontal dashed lines are the correlation lengths £ X1 ( X2 ) 
from Eq. ( 39 1 , which are induced by the finite-entanglement 



cut-off at the critical point. The color shaded background in- 
dicates the two different regimes: Blue is the regime in which 
the iMPS is converged to the exact ground state using bond 
dimensions \i i X2 and red is the scaling regime which shrinks 
with increasing \ ( see main text for further details). 



the critical point. In this regime £ — > £ p hys is independent 
of Xj and all other observables are converged in x- If wc 
measure the critical quantities c and K using Eq. ( 40 ) 
and Eq. (42), we find they renormalize to zero. 



In summary, in the finite entanglement scaling regime, 
1 <C i x <C £phys, we expect Eqs. (40) and p2| to pro- 



duce the critical values c and K, but as £phys <C £ x , the 
MPS converges and we cross over to the true, non-critical 
values c = K = 0. 

The crossover can be analyzed by the following gen- 
eral finite entanglement scaling form. Let g l be a set 
of physical parameters, such as the coupling constants 
or the physical dimension of the system, and let F be a 
scaling observable. Near criticality, F(g) has a scaling 
form determined by the scaling dimensions of F and g % . 
When the system is approximated by an MPS of finite 
X, a new length scale £ x is introduced. While x itself has 
fixed scaling dimension, since £ x ~ x R j we find that it is 
numerically more stable to parameterize the effect of x 
through the MPS correlation length £. We then measure 
F using the MPS, F(g;€). 

The finite entanglement scaling procedure asserts that 
the usual scaling theory still applies to F(g;£), with 
the addition of a single parameter of mass dimension 
[£] = — 1. As usual, the scaling hypothesis allows us 
to rescale F, g, £ in order to eliminate the dependence on 
one parameter. In the usual case in which there are no 
marginal operators, we can linearize the RG equations to 
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determine the scaling dimensions yi, yF, and find 

F(g i ,0 = e ypi F(eyf e g i ,e- e O (43) 

Note that £ is in principle determined both by the 
physical correlation length £ph ys and the finite entangle- 
ment length £ x ; as discussed, at criticality, £ = £ x , while 
at infinite x, £ = £ P hys- Regardless, remains a 

valid coordinate system for the parameters. 



V. DETERMINATION OF THE S = 2 PHASE 
DIAGRAM WITH iDMRG 

In this section we show how the numerical technique 
presented in the previous section was used to investigate 
the Hamiltonian Eq. to establish the phase diagrams 
in Figs. [3] and |4j We begin with a brief summary of the 
characteristics of the phases and the phase transitions in 
that model. We then show representative data used to 
obtain the spin- 2 results of Sec. |III| and discuss in more 
detail its conclusion. 



Characterization of phases and phase 
transitions 



Phases. The phase diagram of the model (IT ) realizes 
several different phases, see for example Figs. 3] and [4] 
For spin-2, the isotropic point, i.e., the Heisenberg model, 
lies in the trivial EH phase in the sense that it does not 
break any symmetry and contains the product state Q. 

The FM and AFM phases are magnetically ordered 
with a nonzero magnetization on each site (S%) 7^ 0. 
The FM phase has nonzero total magnetization (S*) = 
(S* +1 ), while the AFM phase has zero total magnetiza- 
tion (S*) — -(5J +1 ). All the other phases are nonmag- 
netic (S£) = 0. 

The OH phase is an SPTP^HlIlME] stabilized by any of 
the following symmetries: Z 2 x Z 2 rotation symmetry of 
the spins, spatial inversion, and time-reversal symmetry. 
In the presence of Z 2 x2 2 , the phase is characterized by 
a nonlocal string order SO parameter^, 



SO(m,n) = (i> \S^e 



(44) 



which approaches a finite, nonzero value in the thermody- 
namic limit SO = lim| n _ m |_ i . 00 SO(m, n). We also make 
use of another non-local order parameter, Ox, which is 
based on the symmetry under spatial inversion and has 
been introduced in Ref. 1611 This order parameter is 
basically a topological invariant which tells us to which 
cohomology-class the ground state belongs. It is directly 
obtainable from an iMPS in canonical form. The or- 
der parameter Ox takes on values +1/— 1 in the EH/OH 
phases respectively, and thus gives a clear distinction be- 
tween these two phases. If the inversion symmetry is 
broken, the order parameter is set to 0. 



The field theory describing the XY phase is that of 
a free boson <j>. The XY order parameter, which has a 
"stiffness" that determines the energetic cost of modulat- 
ing <f>, is parameterized here by the Luttinger parameter 
K . Within the XY phase K > 0.5, leading to contin- 
uously varying exponents for the correlation functions, 
including one with a power-law decay 



(S m S, N 



i-a(A,D 2 ,D i ) 



(45) 



with correlation length £1 and S„ = S* ± iS%. The ex- 
ponent a(A, D 2 , D4) varies slowly over the XY phase. 
The (S^S*) correlation function is exponentially decay- 
ing, with correlation length £o-^ In some parts of the 
phase diagram, the S z correlation function becomes a 
power-law, while the S ± is exponentially decaying. For 
the parameters we focus on, A > 0, D2 > 0, it is always 
the S correlation function (45) which is relevant. The 



excitation spectrum is gaplcss and the central charge in 
the entire XY phase is c = 1, implying that the entan- 
glement entropy Se diverges throughout the XY phase. 

Phase transitions The AFM O EH phase transition is 
second order, except at large A where it turns first order. 
It is characterized by a continuous vanishing of the mag- 
netic order parameter of the AFM phase and a diverging 
correlation length and entanglement entropy connected 
by c = 1/2. The first order transition is very different 
from the other phase transitions considered here. It does 
not have a critical behavior, no divergence in Se, £ or F, 
and the critical parameters c and K are zero. However, 
it is characterized by a discontinuous slope of the ground 
state energy, from the crossing of two energy levels, and 
a discontinuous jump in the magnetic order parameter. 
The AFM «-» OH phase transition behaves in a similar 
way with the same characteristics along with a vanishing 
string order parameter of the OH phase. 

The EH <H- OH transition is a Gaussian transition de- 
scribed by a conformal field theory with central charge 
c = 1. In addition, the string order parameter (44) van- 



ishes continuously and Ox changes abruptly when enter- 
ing the EH phase. 

The phase transitions discussed so far are between 
gapped phases and the observables obtained from the 
DMRG output normally give a distinct signature at 
the transition. Transitions into the gapless XY phase 
are generally numerically more demanding. Neither the 
EH nor the XY phase has an order parameter with a 
nonzero value in the thermodynamic limit, but transi- 
tions from the OH phase can in principle be determined 
from the vanishing of the string order parameter. In- 
stead, one needs to rely on the characteristics of the tran- 
sitio n, wh ich is of Berezinskii-Kosterlitz-Thouless (BKT) 
£yp e |62i£3] ^ e stiffness weakens, at a critical value of 
K = 0.5 vortices in the </> field unbind, the system be- 
comes disordered and K renormalizes to zero. At the 
BKT transition K = 0.5, however, logarithmic correc- 
tions to scaling are expected due to the presence of a 
marginal operator, which makes the transition difficult 
to study with finite size (or, with iDMRG, finite x) tech- 
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niques. By measuring the stiffness K and the central 
charge c, the BKT point can be determined by a sharp 
drop in K at the value K = 0.5 and a sharp drop in c 
from 1 (XY) to (EH). 



B. Results 



IV 



The results presented below were all obtained with the 
iDMRG algorithm described in Sec 
variety of bond dimensions 50 < x 
specific data set below. The effective correlation length 
obtained at criticality is £ x < 6000. 



< 



We used a wide 
1600, see each 



1. AFM o EH 

An example of data used to obtain the AFM EH 
phase boundary is presented in Fig. [T2| T he local mag- 
netization I (S*) I is obtained from Eq. ( 13 ), the entangle- 
ment entropy Se from Eq. ([6]), the correlation length £ 
from Eq. (fl5|) , and the ground state energy Eqs is defined 



in Sec. |IVj D. The location of the 2 nd order phase tran- 
sition in (a)-(c) shifts with increasing \ and approaches 
the value D^™" EH (A = 1,Z> 4 = 0) = -0.0046 ± 0.0003 
at large x- However, both this shift and the scaling re- 
gion (as defined in Fig. 11) are much smaller than the 
narrow width of the EH phase around the Heisenberg 
point. This phase can hence be detected even with small 
X- Note though, that the shift of the phase transition 
location is large enough to give, at certain values of D 2 , 
an increasing \ (S*)\ with decreasing x, instead of the de- 
creasing trend normally expected in the scaling region. 

At A w 3.8 the transition turns first order. In this case, 
the transition location is obtained from the kink in the 
ground state energy Eqs and the discontinues vanishing 
of the magnetic order parameter, see Fig.[l2|d). All data 
is well converged in x m this type of phase transition 
without critical behavior. 

One of the main questions of this paper, whether the 
Heisenberg point is adiabatically connected to the large- 
ly region, relies on an accurate determination of the nar- 
row part of the EH phase. In addition to locating the 
phase boundaries in Fig. [3j it is important to show that 
the Heisenberg point and the large-D region do not be- 
long to two different phases. To rule out a phase transi- 
tion occurring between them, we note that while the scal- 
ing region surrounding a phase transition can be narrow 
in parameter space, a noticeable increase in quantities 
like Se and £ p hys is normally seen far away from a tran- 
sition. Repeating simulations as in Fig. [l2|a)-(c) across 
the EH phase along lines spaced 8 A = 0.1 apart, we have 
found no such signatures. In particular, the minimum of 
Se and £ p hys remained roughly constant and at the same 
distance away from the AFM phase, all the way up to 
the change in phase transition type at A « 3.8. Further- 
more, no kink in Eqs was observed. With no signs of 




FIG. 12. (Color online) AFM <(->• EH phase transitions, (a)-(c) 
Example of a 2 nd order phase transition. Data as a function 
of D2 away from the Heisenberg point, see Fig. [3] (a) The 
magnetic order parameter |(5^}| of the AFM phase, (b) The 
entanglement entropy Se- (c) The £0 correlation length, (d) 
Example of a 1 st order phase transition. Data along the line 
D 2 = 3.95 - D 4 at A = 4.5, see Fig. [l^a). The magnetic 
order parameter | (5*^)1 (black symbols) and the ground state 
energy _Egs + 4 (red symbols) are plotted. 



a phase transition, an adiabatic connection between the 
Heisenberg point and the large- D region is likely. 



2. EH o XY 

Examples of the BKT type EH <H> XY phase transition, 
located with finite entanglement scaling, are shown in 
Fig. [I3] The central charge c and the stiffness K are 
obtained from a linear fit of the data at various values 
of x to Eqs. (40) and (42) respectively. Examples of 



the Se data can be found in Fig. 12 'b) and the data 
for £ = £1 3> £0 in the XY phase and in the scaling 
region surrounding it behaves in a similar way to £0 in 
Fig. 12 'c) in these regions. The nonzero values of c and 
K in the scaling region fall off continuously from their 
critical values to zero, both with increasing x an d with 
deviation of the coupling constant from criticality. The 
decrease is smallest near criticality and becomes more 
pronounced further away. However, small deviations in 
c and K from their critical values can be observed from 
data simulated with strict convergence criteria. 

Ideally, the phase transition is located by noting that 
in the EH phase c and K scale to zero with increasing x, 
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Method 


d bh-xy 


Df- EH 


20S+iDMRG 


0.045 ± 0.002 


2.42 ±0.05 


LS+ED 


0.043 


2.39 


IOS+DMRG 


0.04 ±0.02 


3.0 ±0.1 



TABLE I. Location of the BKT phase transition as a func- 
tion of D2 away from the Heisenberg point, calculated with 
three different methods. Our results are obtained from two 
observable scaling (20S) and iDMRG, LS±ED results from 
ref. 34 and one observable scaling (IOS) and DMRG results 
from Refs. [31] and 32 



while they scale to their critical values in the XY phase. 
This observation is complicated by the fact that just in- 
side the XY phase the critical parameters, c and espe- 
cially K, have larger values than expected. This devia- 
tion decreases with increasing \ an d disappears further 
into the critical phase. A similar observation was made 
by Song et alP^l for K at a BKT transition in a finite sys- 
tem, and was attributed to the presence of non-universal 
marginal operators at the phase transition causing finite 
size corrections to K. Similar marginal operators could 
give rise to finite entanglement correction in our case, 
however, in our data at the BKT transition these are 
small, due to the large x used in the simulation. Never- 
theless, due to these corrections and the relatively large 
variation in K inside the XY phase, calculation of the 
central charge allows for a more accurate determination 
of the phase boundary than the stiffness calculations. 

The BKT phase transition is more accurately deter- 
mined with scaling in two observable that both depend 
on Xj compared to scaling of a single observable with 
X- Examples of the former include c or K versus £ as 
we use here. Example of the latter are Se or £0, [see 



data in Fig. 12 'b)-(c)], and energy gap SE cx l/£ as used 
in earlier DMRG studies^^Sl w here a finite size scaling 
in L was also required. Another method recently used 
to obtain the phase boundary in finite systems is level- 
spectroscopy (LS), which utilizes the crossing of excited 
energy levels in different magnetization sector s, with 
and without twisted boundary conditions! 34 * 64 * In this 
method, an accurate determination of the BKT phase 
transitions is possible, even using finite systems with ex- 
act diagonalization (ED)P^2£D A comparison between our 
results for the phase transition location in Figs. [13^ a)- (b) 
and those of prior studies is given in Table [T] Good agree- 
ment is observed between the three different studies. 

Outside the XY phase, on the large D 2 side, see Fig.|3j 
the scaling region extends further than at small D 2 , as 
can be seen in Fig. 13 'c) (note the different scale on the 
D 2 axis). A good agreement is again obtained between 
our results and the LS±ED study, see Table [1} The agree- 
ment with the earlier DMRG studies is, however, not as 
good. The large scaling region, makes scaling in one ob- 
servable, like 5E versus Xj more unreliable than scaling 
in two observables and this accounts for most of the dis- 



crepancy. 

Away from the Heisenberg point with increasing A, the 
critical XY phase shrinks in the (D 2 ,D4) planes (see for 
example Fig. EJ). K > 0.5 is only found for A < 2.2, re- 
gardless of the on-site anisotropy strength. In the D4 = 
plane, stiffness calculations have K = 0.5 extending to 
A ps 2.18. This agrees well with the LS±ED study of 
Ref. [321 An example of the decrease in K when leav- 
ing the XY phase along the line D 2 = 0.85A — 0.055 at 
D4 = can be seen in Fig. 13 'd). Apart from the slow 
decrease of K across the K = 0.5 line, no other clear 
features can be observed, maybe with the exception of 
a finite x scaling around K — 0.5. However, this scal- 
ing is smaller than at the previous two BKT transitions 
we considered in Fig. [13] (data for K not shown for the 
second of those), indicating that the XY phase could be 
slightly smaller in simulations with larger x m this re- 
gion. No decrease in K with x can be observed outside 
the XY phase. The reason is that our data can not dis- 
tinguish these states from critical states with c = 1, due 
to the nearby (in D4) EH f-> OH phase transition. We 
discuss this in more detail below. Earlier DMRG results 
from the vanishing of the energy gap had a much larger 
XY phase, extending around this line, to A ps 3.8221. xhe 
close presence of the EH-OH phase transition is also the 
reason for the much larger XY phase that was obtained 
in earlier DMRG calculations. 



3. EH <rt OH 

The second main question for the spin-2 chain we ad- 
dress, is whether the OH phase exists in the model ([I]). 
We argued in Sec. [TT] that with the addition of the D4 
term one effectively projects the system to a spin-1 sys- 
tem, realizing the corresponding spin-1 phase diagram. 
Indeed, as shown in Fig. HJ with increasing D4 an OH 
phase with ground states that can be described accu- 
rately with an iMPS is observed in a large part of the 
phase diagram. The ground states have doubly degen- 
erate entanglement spectra, non-zero string order, and 
Ox = — 1. Below, we focus on exploring how close to 
the XXZ-chain the OH phase extends. Hence, we will 
primarily investigate the EH -h- OH phase transition for 
small on-site anisotropics, which is the experimentally 
most relevant case. 

The OH phase extends to A OH " EH ps 5.3 close to the 
AFM phase, much further than the XY phase which only 
extends toA« 2.2. For example, the relevant part of the 
(Z?2,£>4) phase plane at A 



4.5 is shown in Fig. 14 a). 
The OH phase extends almost all the way to the D4 
plane. In fact, the transition into the EH phase at 
£,eh-oh = 0.0135 ± 0.0010 for A = 4.5 and D 2 = 3.95 is 
approximately the closest approach of the OH phase to 
the D4 = plane for any parameters in Eq. (JlJ. 

The location of the EH f-> OH phase transition can be 
accurately determined from many different observables 



far away from the critical phase. Fig. 14 'b) shows, for 
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FIG. 13. 

obtained 
X values 
(a) were 
D2 away 
show the 
Stiffness 
Fig. E 



(Color online) BKT phase transition. The data are 
from finite entanglement scaling averaged over the 
given in the insets [for (b) the same x values as in 
used], (a)-(c) Critical quantities as a function of 
from the Heisenberg point, see Fig. [3] (a) and (c) 
central charge, while (b) and (d) the stiffness, (d) 
along the line D 2 = 0.85 • A - 0.055 at D 4 = 0, see 



example, how the string order [obtained from Eqs. ( 13 ) 
and ( 14 ) in the limit x ~^ 00 ] and the central charge vary 
across the phase transition along the line D 2 = 3.95 — D 4 
[see black dashed line in Fig. [Tl^a)]. The string order 



has a small value in the OH phase, and abruptly changes 
to zero at the phase transition. A peak in c to a value 
slightly larger than 1, is observed at the phase transition 
consistent with the expected critical behavior. Since in 
this data in the OH phase, we are never very far from 
the phase transition, the decrease of c with \ is slow. 
The OH phase disappears again from this phase diagram 
[Fig. fl4{a)] at large D 4 in the D 4 — > 00 limit. 

Closer to the critical XY phase, at smaller values 
of A, it is harder to accurately locate the EH 4-» OH 
phase transition. The peaks in the critical properties 
that diverge in the thermodynamic limit get broader and 
broader, for the x we can simulate, while moving to 
a slightly larger D4 as A is decreased. For A < 2.6 
it is unclear where the peak is and close to the XY 



phase 



1 for -0.03 < £> 4 < 0.12. The order pa- 



rameters of the OH phase (SO and projected inversion 
symmetry) locate the phase transition more accurately. 
Fig. [l4^c) shows Ox, calculated as in Ref. [6U along the 
line D 2 — 2.155 — D 4 at A = 2.6. The phase transition is 
obtained at Df H -° H (A = 2.6, D 2 = 2.105) = 0.05 ±0.04. 



No scaling with x is observed for this order parameter. 

Even closer to the XY phase both the expectation value 
and the uncertainty of the location of the EH ■<-» OH 
phase transition increase, the former to larger D 4l from 
similar data sets as in Fig. [l4^c) (not shown). However, 
the uncertainty is always smaller than the distance to the 
D 4 = plane as can be seen in Fig. [l4|d), where no sign 
of an OH phase can be seen. The data is along the line 
D 2 = 0.85 • A — 0.055 at D 4 = which is in the center of 
where Refs. 35 and 36 found a narrow OH phase. Also, 



the stiffness in Fig. 13 d) was calculated along this line. 



Note, that no sign of the XY O EH phase transition 
can be seen in the Ox order parameter. As the critical 
state cannot be represented exactly with a finite dimen- 
sional MPS, the order parameter Ox is not well defined 
here. The DMRG optimization of the finite dimensional 
MPS, however, yields a state which has an approximate 
inversion symmetric ground state which lies in the EH 
phase with Ox = 1- 

Parts of the phase diagram of the spin-2 XXZ-chain 
with on-site anisotropy are difficult to study numerically. 
This is mainly due to the presence of the critical XY 
phase, the EH «-» OH phase transition which has some 
of the same critical properties (c = 1), and the large 
scaling region that surrounds them. This is especially 
true around the line where these two phase transition 
types meet. Earlier DMRG studieP^t^H could not de- 
termine the location of the phase boundaries accurately 
where the scaling region is large. While increased com- 
putational power helps, scaling in two observables that 
both depend on \ rather than a direct scaling with x is 
even more important. In all other parts of the phase di- 
agram our results are in good agreement with the above 
mentioned studies. 

The location of the BKT transition agrees well with 
that obtained in the LS+ED studies^^U. S ome d iscrep- 
ancy between our result and LS the results^^fil can hg 
found in the location of the EH <H- OH phase transition 
close to the XY phase, but the difference in the obtained 
values of D 4 is nevertheless small. The EE 35 study con- 
sidered small systems (L < 12), and obtained a roughly 
linear scaling in 1/L 2 , just as what was obtained for 
the BKT phase transition. However, for larger systems 
(L < 28)P a different scaling was observed in DMRG 
leading to a smaller OH phase. In the region with three 
nearby competing phases and a huge correlation length 
£phys ^ 10000, scaling from small systems can sometimes 
be unreliable. Even some of our data is uncertain in this 
region, since a bond dimension of x < 560 (800 for some 
points) is not enough to clearly distinguish between the 
many nearly degenerate energy eigenstates present. 

Finally, we mention that a recent study^ on a related 
model showed that the OH phase also appears in a spin-2 
chain at the SO(5) symmetric point, obtained by tuning 
the J p 's in H = £„ £j=i J p (S n ■ S n+1 )P. 
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FIG. 14. (Color online) (a) Phase diagram at A = 4.5. (b) 
String order SO and central charge c across the EH-OH phase 
transition along the line Di = 3.95 — Di shown in (a), (c) 
Projected inversion symmetry along the line D2 = 2.155 — D4 
at A = 2.6. (d) Projected inversion symmetry Ox along the 
line D 2 = 0.85 • A - 0.055 at D A = 0, see Fig. [3] 



VI. SPIN-3 XXZ-CHAIN 

Having applied the iDMRG algorithm to study the 
phase diagram of the spin-2 XXZ Heisenberg model, we 
provide the first, to our knowledge, numerical investiga- 
tion of the spin-3 XXZ-chain ^ around the Heisenberg 
point. The approach to the classical limit S 00 is 
analyzed by comparing the results with the 5 = 1 and 2 
cases. The spin-3 chain is hard to investigate numerically, 
not primarily because of the large local Hilbert space, but 
due to the tiny OH phase with the correspondingly large 
correlation lengths, see Sec. |TTj 

In the semi-classical limit of large but finite S, an ex- 
ponentially decreasing gap is obtained at the Heisenberg 
point: SEg ~ S 2 e"^ s as S — > 00, with a spin wave 
velocity v = 6Es£ = 2S relating it to the correlation 
lengthPZI The width of the Haldane phase in A also de- 
crease rapidly with S in the semi-classical limit. 

For 5" = 1,2 the size of the AFM phase is slightly 
overestimated at finite x but decreases with increasing 
X, see Fig. 12 'a). The same occurs for S — 3, see 
Fig. 15 [a), which plots the value A* of A at which the 
AFM magnetic order parameter vanishes ((S^) — > 0) 
as a function of 1/%, for \ < 1120. In the ther- 



modynamic limit, the OH «-» AFM phase transition is 
at AgS^FM = 1.000045 ± 0.000020, very close to the 
Heisenberg point. The shift of the phase boundary with 
increasing x is also observed in the peak of the £0 corre- 



Spin 


^XY-OH 


^OH-AFM 


1 


0.0 


1.19 ±0.01 


2 


0.962 ± 0.004 


1.0039 ± 0.0005 


3 


0.99965 ± 0.00010 


1.000045 ± 0.000020 



TABLE II. Location of the phase transitions in and out of the 
Haldane phase in S = 1,2,3 XXZ-chains. For S = 1,2 this 
data has been obtained before, for example in Ref. 1401 They 
are given here within a wide interval for a comparison to the 
spin-3 results. 



aXY-OH 



lation length, see Fig. 15 Id). 

The XY <-> OH phase transition at 
0.99965±0. 00010 is located by the vanishing of the string 
order (not shown) and central charge calculations, see 
Fig. [l5{c) . Locating the point at which the string order 
vanishes is challenging. The string order decays with a 
power-law in the critical phase and its value in the ther- 
modynamic limit in the gapped OH phase approaches 
zero exponentially at the XY «-» OH phase transition. 
However, with a careful scaling of the string order versus 
X the same phase transition location as with the central 
charge is obtained. The projected inversion symmetry 
Ox is not useful here since it is undefined at critical- 
ity and can not detect this type of phase transition, see 
Fig. 14 'd) and its discussion. 

At the Heisenberg point, the OH phase is obtained as 
the ground state for x ^ 800, as can be seen from the 
point at which the curve in Fig. [l5|^a) crosses A* = 1. 
The correlation length £ s = 3 = 520 ± 60 and the string 



order SO s " 3 = 0.162 ± 0.002, see Fig. [llgd), in the 
thermodynamic limit is hence obtained from the data 
at x = 800, 1120, 1340, 1600 by scaling. Comparing with 
S = 1,2 where £ 5=1 = 5.77±0.01 and £ s = 2 = 47.2 ±0.5, 
we indeed have an exponential growth, but with a bit 
smaller exponent than in the asymptotic limit S — > 00. 
Nevertheless, it is large enough to make the Haldane 
phases in the quantum S > 4 XXZ-chains unreachable 
with the current approach. Moreover, the width of the 
Haldane phase decrease even faster than l/£, with the 
XY o Haldane phase transition A^ Y " H and the Haldane 
«-> AFM phase transition Af " AFM presented in Table [n] 
for S — 1,2,3. For spin-3 the narrow OH phase is thus 
not observed for x 200. 



VII. CONCLUSIONS 

In this work we have carefully studied the phase dia- 
gram of the spin-2 XXZ Heisenberg model with on-site 
anisotropies using the iDMRG algorithm. We have es- 
tablished that the SPTP OH phase is obtained for not 
too large values of the on-site anisotropy, and that the 
Heisenberg point and the large-D region are adiabati- 
cally connected and belong to the EH phase. We have 
also provided the first, to our knowledge, study of the 
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corresponding spin-3 XXZ Heisenberg model. Thereby 
we have numerically further verified the exponentially 
decreasing gap and size of the Haldane phase with in- 
creasing spin. Many of these conclusions were greatly 
aided by scaling in two observables that both vary with 
the bond dimension x, which can locate deviations from 
criticality more accurately than scaling directly with x- 
Additionally, we have provided a self-contained and di- 
dactic introduction to the iDMRG algorithm, which was 
applied in this study. 
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